setwd("[your home directory]/IRF_Figures/InputData")
rm(list = ls())
library(tidyselect)
library(readr)
library(lpirfs)
library(gridExtra)
library(ggpubr)
library(cowplot)
library(colorspace)
library(foreign)
library(haven)

#-----------------------------Clean Data For Charts-------------------------####
mindate=1
maxdate=92 # 121=only run through 2019, 92=run thru 2012q2 

dfa_gini=read_csv("dfa_gini_networth_low.csv")
dfa_networth_shares <- read_csv("dfa_shares_bus_cycle2_networth_2019q3.csv") 

NW_shares=dfa_networth_shares['networth']

Group_shares=t(matrix(unlist(NW_shares), nrow=6))
Group_shares=Group_shares[mindate:maxdate,]

TOP01=data.frame(Group_shares[,1:1])
#TOP1=data.frame(rowSums(Group_shares[,1:2]))  #combines top .1% and top 1%
#TOP10=data.frame(rowSums(Group_shares[,1:3])) #combines top .1%, top 1% and top 10%
TOP1=data.frame(Group_shares[,2:2])
TOP10=data.frame(Group_shares[,3:3])
BOT90=data.frame(rowSums(Group_shares[,4:6]))
P7090=data.frame(Group_shares[,4:4])
P5070=data.frame(Group_shares[,5:5])
BOT50=data.frame(Group_shares[,6:6])

DFAGINI=dfa_gini['gini']
GINI=DFAGINI[mindate:maxdate,]

FFrate <- read.csv("FEDFUNDS.csv")
FFrate=FFrate[mindate:maxdate,]
FFR=FFrate['FEDFUNDS']

# upload_shocks_june2023.csv -- includes the Gertler and Karadi 2015 shocks
MPShocks_all<- read.csv("upload_shocks_june2023.csv")
MPS_RR_shocks=MPShocks_all['RR_Original']
MPS_RR_shocks_ext=MPShocks_all['RR_Extended']
MPS_NS_shocks1=MPShocks_all['NS_FFR']
MPS_NS_shocks2=MPShocks_all['NS_News']
MPS_BRW_shocks=MPShocks_all['BRW']

# Romer and Romer
MPS_RR=as.data.frame(MPS_RR_shocks[mindate:maxdate,])
MPS_RR_ext=as.data.frame(MPS_RR_shocks_ext[mindate:maxdate,])
# Nakamura Steinsson
MPS_NS1=as.data.frame(MPS_NS_shocks1[mindate:maxdate,])
MPS_NS2=as.data.frame(MPS_NS_shocks2[mindate:maxdate,])
MPS_BRW=as.data.frame(MPS_BRW_shocks[mindate:maxdate,])

# New Gertler shocks (June 2023)
#MPShocks_Gert<- read.csv("factor_data_qtrly.csv")
MPS_GK_shock1=MPShocks_all['GK_MP1']
MPS_GK_shock2=MPShocks_all['GK_FF4']

MPS_GK1=as.data.frame(MPS_GK_shock1[mindate:maxdate,])
MPS_GK2=as.data.frame(MPS_GK_shock2[mindate:maxdate,])

gdp<- read.csv("GDPGAP.csv")
gdp=gdp[mindate:maxdate,]
GDP=gdp['GDPC1_GDPPOT']

ugap<- read.csv("ugap.csv")
ugap=ugap[mindate:maxdate,]
UGAP=ugap['NROU_UNRATE']

sp500 <- read.csv("upload_sp500.csv")
sp500 =sp500['Return_ann']

sp=sp500[(mindate:maxdate),]
SP500=sp

STOCK= sp


### GK1 Shock GINI DONE#####################################
shock="MP_GK1"
IRF_data_gini_exog_data=as.data.frame(cbind(STOCK,  GDP, UGAP))
IRF_data_gini_endog_data=as.data.frame(cbind(GINI, -FFR))

results_lin_iv_gini_GK1 <- lp_lin_iv(endog_data=IRF_data_gini_endog_data,
                                     lags_endog_lin = 4,
                                     shock     = -MPS_GK1,
                                     exog_data = IRF_data_gini_exog_data,
                                     lags_exog =2,
                                     trend          = 0,
                                     confint        = 1.645,
                                     hor            = 16)



### NS1 Shock GINI DONE####
shock="MP_NS1"
IRF_data_gini_exog_data=as.data.frame(cbind(STOCK,  GDP, UGAP))
IRF_data_gini_endog_data=as.data.frame(cbind(GINI, -FFR))

results_lin_iv_gini_NS1 <- lp_lin_iv(endog_data=IRF_data_gini_endog_data,
                                     lags_endog_lin = 4,
                                     shock     = -MPS_NS1,
                                     exog_data = IRF_data_gini_exog_data,
                                     lags_exog =2,
                                     trend          = 0,
                                     confint        = 1.645,
                                     hor            = 16)


### GK1 Shock Top .1 percent DONE#####################################
shock="MP_GK1"
IRF_data_top01_exog_data=as.data.frame(cbind(STOCK,  GDP, UGAP))
IRF_data_top01_endog_data=as.data.frame(cbind(TOP01, -FFR))

results_lin_iv_top01 <- lp_lin_iv(endog_data=IRF_data_top01_endog_data,
                                  lags_endog_lin = 4,
                                  shock     = -MPS_GK1,
                                  exog_data = IRF_data_top01_exog_data,
                                  lags_exog =2,
                                  trend          = 0,
                                  confint        = 1.645,
                                  hor            = 16)


### GK1 Shock Top 1 percent DONE####
shock="MP_GK1"
IRF_data_top1_exog_data=as.data.frame(cbind(STOCK,  GDP, UGAP))
IRF_data_top1_endog_data=as.data.frame(cbind(TOP1, -FFR))

results_lin_iv_top1 <- lp_lin_iv(endog_data=IRF_data_top1_endog_data,
                                 lags_endog_lin = 4,
                                 shock     = -MPS_GK1,
                                 exog_data = IRF_data_top1_exog_data,
                                 lags_exog =2,
                                 trend          = 0,
                                 confint        = 1.645,
                                 hor            = 16)


### GK1 Shock Top 10 percent DONE####
shock="MP_GK1"
IRF_data_top10_exog_data=as.data.frame(cbind(STOCK,  GDP, UGAP))
IRF_data_top10_endog_data=as.data.frame(cbind(TOP10, -FFR))

results_lin_iv_top10 <- lp_lin_iv(endog_data=IRF_data_top10_endog_data,
                                  lags_endog_lin = 4,
                                  shock     = -MPS_GK1,
                                  exog_data = IRF_data_top10_exog_data,
                                  lags_exog =2,
                                  trend          = 0,
                                  confint        = 1.645,
                                  hor            = 16)


### GK1 Shock 70 to 90 percentile DONE####
shock="MP_GK1"
IRF_data_p7090_exog_data=as.data.frame(cbind(STOCK,  GDP, UGAP))
IRF_data_p7090_endog_data=as.data.frame(cbind(P7090, -FFR))

results_lin_iv_p7090 <- lp_lin_iv(endog_data=IRF_data_p7090_endog_data,
                                  lags_endog_lin = 4,
                                  shock     = -MPS_GK1,
                                  exog_data = IRF_data_p7090_exog_data,
                                  lags_exog =2,
                                  trend          = 0,
                                  confint        = 1.645,
                                  hor            = 16)


### GK1 Shock 50 to 70 percentile DONE####
shock="MP_GK1"
IRF_data_p5070_exog_data=as.data.frame(cbind(STOCK,  GDP, UGAP))
IRF_data_p5070_endog_data=as.data.frame(cbind(P5070, -FFR))

results_lin_iv_p5070 <- lp_lin_iv(endog_data=IRF_data_p5070_endog_data,
                                  lags_endog_lin = 4,
                                  shock     = -MPS_GK1,
                                  exog_data = IRF_data_p5070_exog_data,
                                  lags_exog =2,
                                  trend          = 0,
                                  confint        = 1.645,
                                  hor            = 16)


### GK1 Shock Bottom 50 percentile DONE####
shock="MP_GK1"
IRF_data_bot50_exog_data=as.data.frame(cbind(STOCK,  GDP, UGAP))
IRF_data_bot50_endog_data=as.data.frame(cbind(BOT50, -FFR))

results_lin_iv_bot50 <- lp_lin_iv(endog_data=IRF_data_bot50_endog_data,
                                  lags_endog_lin = 4,
                                  shock     = -MPS_GK1,
                                  exog_data = IRF_data_bot50_exog_data,
                                  lags_exog =2,
                                  trend          = 0,
                                  confint        = 1.645,
                                  hor            = 16)


#--------------------------------Figure 13 Panel A--------------------------####
figure_13a <- function()
{
  linear_plots_gini_exog_MPGK1<- plot_lin(results_lin_iv_gini_GK1)
  linear_plots_gini_exog_MPGK1[[1]][["coordinates"]][["limits"]][["y"]]=c(-.025,.05)
  linear_plots_gini_exog_MPGK1[[1]][["labels"]][["title"]]="GK1 / GINI"
  linear_plots_gini_exog_MPGK1[[1]][["labels"]][["y"]]="GK1 Effect on Wealth Gini Coefficient"
  linear_plots_gini_exog_MPGK1[[1]][["labels"]][["x"]]="Quarters"
  linear_plots_gini_exog_MPGK1[[1]][["theme"]][["axis.title.y"]][["size"]]=c(10,10)
  linear_plots_gini_exog_MPGK1[[1]][["theme"]][["axis.title.x"]][["size"]]=c(10,10)
  linear_plots_gini_exog_MPGK1[[1]][["theme"]][["axis.text"]][["colour"]]='black'
  linear_plots_gini_exog_MPGK1[[1]][["theme"]][["axis.text"]][["size"]]=10
  linear_plots_gini_exog_MPGK1[1]
}
figure_13a()

#--------------------------------Figure 13 Panel B--------------------------####
figure_13b <- function()
{
  linear_plots_gini_exog_MPNS1<- plot_lin(results_lin_iv_gini_NS1)
  linear_plots_gini_exog_MPNS1[[1]][["coordinates"]][["limits"]][["y"]]=c(-.010,.015)
  linear_plots_gini_exog_MPNS1[[1]][["labels"]][["title"]]="NS1 / GINI"
  linear_plots_gini_exog_MPNS1[[1]][["labels"]][["y"]]="NS1 Effect on Wealth Gini Coefficient"
  linear_plots_gini_exog_MPNS1[[1]][["labels"]][["x"]]="Quarters"
  linear_plots_gini_exog_MPNS1[[1]][["theme"]][["axis.title.y"]][["size"]]=c(10,10)
  linear_plots_gini_exog_MPNS1[[1]][["theme"]][["axis.title.x"]][["size"]]=c(10,10)
  linear_plots_gini_exog_MPNS1[[1]][["theme"]][["axis.text"]][["colour"]]='black'
  linear_plots_gini_exog_MPNS1[[1]][["theme"]][["axis.text"]][["size"]]=10
  linear_plots_gini_exog_MPNS1[1]
}
figure_13b()

#--------------------------------Figure 14 Panel A--------------------------####
figure_14a <- function()
{
  linear_plots_top01_exog_MPGK1<- plot_lin(results_lin_iv_top01)
  linear_plots_top01_exog_MPGK1[[1]][["coordinates"]][["limits"]][["y"]]=c(-.025,.1)
  linear_plots_top01_exog_MPGK1[[1]][["labels"]][["title"]]="GK1 / Top .1%"
  linear_plots_top01_exog_MPGK1[[1]][["labels"]][["y"]]="GK1 Effect on 0.1%"
  linear_plots_top01_exog_MPGK1[[1]][["labels"]][["x"]]="Quarters"
  linear_plots_top01_exog_MPGK1[[1]][["theme"]][["axis.title.y"]][["size"]]=c(10,10)
  linear_plots_top01_exog_MPGK1[[1]][["theme"]][["axis.title.x"]][["size"]]=c(10,10)
  linear_plots_top01_exog_MPGK1[[1]][["theme"]][["axis.text"]][["colour"]]='black'
  linear_plots_top01_exog_MPGK1[[1]][["theme"]][["axis.text"]][["size"]]=10
  linear_plots_top01_exog_MPGK1[1]
}
figure_14a()

#--------------------------------Figure 14 Panel B--------------------------####
figure_14b <- function()
{
  linear_plots_top1_exog_MPGK1<- plot_lin(results_lin_iv_top1)
  linear_plots_top1_exog_MPGK1[[1]][["coordinates"]][["limits"]][["y"]]=c(-.08,.05)
  linear_plots_top1_exog_MPGK1[[1]][["labels"]][["title"]]="GK1 / Top 99-99.9%"
  linear_plots_top1_exog_MPGK1[[1]][["labels"]][["y"]]="GK1 Effect on Top 99-99.9%"
  linear_plots_top1_exog_MPGK1[[1]][["labels"]][["x"]]="Quarters"
  linear_plots_top1_exog_MPGK1[[1]][["theme"]][["axis.title.y"]][["size"]]=c(10,10)
  linear_plots_top1_exog_MPGK1[[1]][["theme"]][["axis.title.x"]][["size"]]=c(10,10)
  linear_plots_top1_exog_MPGK1[[1]][["theme"]][["axis.text"]][["colour"]]='black'
  linear_plots_top1_exog_MPGK1[[1]][["theme"]][["axis.text"]][["size"]]=10
  linear_plots_top1_exog_MPGK1[1]
  
}
figure_14b()

#--------------------------------Figure 14 Panel C--------------------------####
figure_14c <- function()
{
  linear_plots_top10_exog_MPGK1<- plot_lin(results_lin_iv_top10)
  linear_plots_top10_exog_MPGK1[[1]][["coordinates"]][["limits"]][["y"]]=c(-.08,.050)
  linear_plots_top10_exog_MPGK1[[1]][["labels"]][["title"]]="GK1 / Top 90-99%"
  linear_plots_top10_exog_MPGK1[[1]][["labels"]][["y"]]="GK1 Effect on Top 90-99%"
  linear_plots_top10_exog_MPGK1[[1]][["labels"]][["x"]]="Quarters"
  linear_plots_top10_exog_MPGK1[[1]][["theme"]][["axis.title.y"]][["size"]]=c(10,10)
  linear_plots_top10_exog_MPGK1[[1]][["theme"]][["axis.title.x"]][["size"]]=c(10,10)
  linear_plots_top10_exog_MPGK1[[1]][["theme"]][["axis.text"]][["colour"]]='black'
  linear_plots_top10_exog_MPGK1[[1]][["theme"]][["axis.text"]][["size"]]=10
  linear_plots_top10_exog_MPGK1[1]
}
figure_14c()

#--------------------------------Figure 14 Panel D--------------------------####
figure_14d <- function()
{
  linear_plots_p7090_exog_MPGK1<- plot_lin(results_lin_iv_p7090)
  linear_plots_p7090_exog_MPGK1[[1]][["coordinates"]][["limits"]][["y"]]=c(-.12,.03)
  linear_plots_p7090_exog_MPGK1[[1]][["labels"]][["title"]]="GK1 / Top 70-90%"
  linear_plots_p7090_exog_MPGK1[[1]][["labels"]][["y"]]="GK1 Effect on Top 70-90%"
  linear_plots_p7090_exog_MPGK1[[1]][["labels"]][["x"]]="Quarters"
  linear_plots_p7090_exog_MPGK1[[1]][["theme"]][["axis.title.y"]][["size"]]=c(10,10)
  linear_plots_p7090_exog_MPGK1[[1]][["theme"]][["axis.title.x"]][["size"]]=c(10,10)
  linear_plots_p7090_exog_MPGK1[[1]][["theme"]][["axis.text"]][["colour"]]='black'
  linear_plots_p7090_exog_MPGK1[[1]][["theme"]][["axis.text"]][["size"]]=10
  linear_plots_p7090_exog_MPGK1[1]
  
}
figure_14d()

#--------------------------------Figure 14 Panel E--------------------------####
figure_14e <- function()
{
  linear_plots_p5070_exog_MPGK1<- plot_lin(results_lin_iv_p5070)
  linear_plots_p5070_exog_MPGK1[[1]][["coordinates"]][["limits"]][["y"]]=c(-.04,.03)
  linear_plots_p5070_exog_MPGK1[[1]][["labels"]][["title"]]="GK1 / Top 50-70%"
  linear_plots_p5070_exog_MPGK1[[1]][["labels"]][["y"]]="GK1 Effect on Top 50-70%"
  linear_plots_p5070_exog_MPGK1[[1]][["labels"]][["x"]]="Quarters"
  linear_plots_p5070_exog_MPGK1[[1]][["theme"]][["axis.title.y"]][["size"]]=c(10,10)
  linear_plots_p5070_exog_MPGK1[[1]][["theme"]][["axis.title.x"]][["size"]]=c(10,10)
  linear_plots_p5070_exog_MPGK1[[1]][["theme"]][["axis.text"]][["colour"]]='black'
  linear_plots_p5070_exog_MPGK1[[1]][["theme"]][["axis.text"]][["size"]]=10
  linear_plots_p5070_exog_MPGK1[1]
}
figure_14e()

#--------------------------------Figure 14 Panel F--------------------------####
figure_14f <- function()
{
  linear_plots_bot50_exog_MPGK1<- plot_lin(results_lin_iv_bot50)
  linear_plots_bot50_exog_MPGK1[[1]][["coordinates"]][["limits"]][["y"]]=c(-.02,.03)
  linear_plots_bot50_exog_MPGK1[[1]][["labels"]][["title"]]="GK1 / Bottom 50%"
  linear_plots_bot50_exog_MPGK1[[1]][["labels"]][["y"]]="GK1 Effect on Bottom 50%"
  linear_plots_bot50_exog_MPGK1[[1]][["labels"]][["x"]]="Quarters"
  linear_plots_bot50_exog_MPGK1[[1]][["theme"]][["axis.title.y"]][["size"]]=c(10,10)
  linear_plots_bot50_exog_MPGK1[[1]][["theme"]][["axis.title.x"]][["size"]]=c(10,10)
  linear_plots_bot50_exog_MPGK1[[1]][["theme"]][["axis.text"]][["colour"]]='black'
  linear_plots_bot50_exog_MPGK1[[1]][["theme"]][["axis.text"]][["size"]]=10
  linear_plots_bot50_exog_MPGK1[1]
}
figure_14f()

#----------------------------------Export Series----------------------------####
setwd("[your home directory]/IRF_Figures/Figures13_14/")
export_figure13_series <- function()
{
  # Create Figure 13a Series 
  figure13aseries_high <- as.data.frame(t(results_lin_iv_gini_GK1$irf_lin_up[, ]))
  figure13aseries_high <- as.data.frame(figure13aseries_high$V1)
  
  figure13aseries_mean <- as.data.frame(t(results_lin_iv_gini_GK1$irf_lin_mean[, ]))
  figure13aseries_mean <- as.data.frame(figure13aseries_mean$V1)
  
  figure13aseries_low <- as.data.frame(t(results_lin_iv_gini_GK1$irf_lin_low[, ]))
  figure13aseries_low <- as.data.frame(figure13aseries_low$V1)
  
  # Create Figure 13b Series 
  figure13bseries_high <- as.data.frame(t(results_lin_iv_gini_NS1$irf_lin_up[, ]))
  figure13bseries_high <- as.data.frame(figure13bseries_high$V1)
  
  figure13bseries_mean <- as.data.frame(t(results_lin_iv_gini_NS1$irf_lin_mean[, ]))
  figure13bseries_mean <- as.data.frame(figure13bseries_mean$V1)
  
  figure13bseries_low <- as.data.frame(t(results_lin_iv_gini_NS1$irf_lin_low[, ]))
  figure13bseries_low <- as.data.frame(figure13bseries_low$V1)
  
  # Merge all series together
  figure13series <- list(figure13aseries_high, figure13aseries_mean, figure13aseries_low, figure13bseries_high, figure13bseries_mean, figure13bseries_low)
  figure13series <- as.data.frame(figure13series)
  figure13series <- (mutate(figure13series, quarter=row_number()))
  quarter <- as.data.frame(figure13series$quarter)
  figure13series <- list(quarter, figure13aseries_high, figure13aseries_mean, figure13aseries_low, figure13bseries_high, figure13bseries_mean, figure13bseries_low)
  figure13series <- as.data.frame(figure13series)
  figure13series <- rename(figure13series, quarters = `figure13series.quarter`, figure13ahigh = `figure13aseries_high.V1`, figure13bhigh = `figure13bseries_high.V1`, figure13amean = `figure13aseries_mean.V1`, figure13bmean = `figure13bseries_mean.V1`, figure13alow = `figure13aseries_low.V1`, figure13blow = `figure13bseries_low.V1`)
  
  # Export the series 
  write.csv(figure13series, "figure13series.csv")
}
export_figure13_series()

export_figure14_series <- function()
{
  # Create Figure 14a Series 
  figure14aseries_high <- as.data.frame(t(results_lin_iv_top01$irf_lin_up[, ]))
  figure14aseries_high <- as.data.frame(figure14aseries_high$V1)
  
  figure14aseries_mean <- as.data.frame(t(results_lin_iv_top01$irf_lin_mean[, ]))
  figure14aseries_mean <- as.data.frame(figure14aseries_mean$V1)
  
  figure14aseries_low <- as.data.frame(t(results_lin_iv_top01$irf_lin_low[, ]))
  figure14aseries_low <- as.data.frame(figure14aseries_low$V1)
  
  # Create Figure 14b Series 
  figure14bseries_high <- as.data.frame(t(results_lin_iv_top1$irf_lin_up[, ]))
  figure14bseries_high <- as.data.frame(figure14bseries_high$V1)
  
  figure14bseries_mean <- as.data.frame(t(results_lin_iv_top1$irf_lin_mean[, ]))
  figure14bseries_mean <- as.data.frame(figure14bseries_mean$V1)
  
  figure14bseries_low <- as.data.frame(t(results_lin_iv_top1$irf_lin_low[, ]))
  figure14bseries_low <- as.data.frame(figure14bseries_low$V1)
  
  # Create Figure 14c Series 
  figure14cseries_high <- as.data.frame(t(results_lin_iv_top10$irf_lin_up[, ]))
  figure14cseries_high <- as.data.frame(figure14cseries_high$V1)
  
  figure14cseries_mean <- as.data.frame(t(results_lin_iv_top10$irf_lin_mean[, ]))
  figure14cseries_mean <- as.data.frame(figure14cseries_mean$V1)
  
  figure14cseries_low <- as.data.frame(t(results_lin_iv_top10$irf_lin_low[, ]))
  figure14cseries_low <- as.data.frame(figure14cseries_low$V1)
  
  # Create Figure 14d Series 
  figure14dseries_high <- as.data.frame(t(results_lin_iv_p7090$irf_lin_up[, ]))
  figure14dseries_high <- as.data.frame(figure14dseries_high$V1)
  
  figure14dseries_mean <- as.data.frame(t(results_lin_iv_p7090$irf_lin_mean[, ]))
  figure14dseries_mean <- as.data.frame(figure14dseries_mean$V1)
  
  figure14dseries_low <- as.data.frame(t(results_lin_iv_p7090$irf_lin_low[, ]))
  figure14dseries_low <- as.data.frame(figure14dseries_low$V1)
  
  # Create Figure 14e Series 
  figure14eseries_high <- as.data.frame(t(results_lin_iv_p5070$irf_lin_up[, ]))
  figure14eseries_high <- as.data.frame(figure14eseries_high$V1)
  
  figure14eseries_mean <- as.data.frame(t(results_lin_iv_p5070$irf_lin_mean[, ]))
  figure14eseries_mean <- as.data.frame(figure14eseries_mean$V1)
  
  figure14eseries_low <- as.data.frame(t(results_lin_iv_p5070$irf_lin_low[, ]))
  figure14eseries_low <- as.data.frame(figure14eseries_low$V1)
  
  # Create Figure 12f Series 
  figure14fseries_high <- as.data.frame(t(results_lin_iv_bot50$irf_lin_up[, ]))
  figure14fseries_high <- as.data.frame(figure14fseries_high$V1)
  
  figure14fseries_mean <- as.data.frame(t(results_lin_iv_bot50$irf_lin_mean[, ]))
  figure14fseries_mean <- as.data.frame(figure14fseries_mean$V1)
  
  figure14fseries_low <- as.data.frame(t(results_lin_iv_bot50$irf_lin_low[, ]))
  figure14fseries_low <- as.data.frame(figure14fseries_low$V1)
  
  # Merge all series together
  figure14series <- list(figure14aseries_high, figure14aseries_mean, figure14aseries_low, figure14bseries_high, figure14bseries_mean, figure14bseries_low, figure14cseries_high, figure14cseries_mean, figure14cseries_low, figure14dseries_high, figure14dseries_mean, figure14dseries_low, figure14eseries_high, figure14eseries_mean, figure14eseries_low, figure14fseries_high, figure14fseries_mean, figure14fseries_low)
  figure14series <- as.data.frame(figure14series)
  figure14series <- (mutate(figure14series, quarter=row_number()))
  quarter <- as.data.frame(figure14series$quarter)
  figure14series <- list(quarter, figure14aseries_high, figure14aseries_mean, figure14aseries_low, figure14bseries_high, figure14bseries_mean, figure14bseries_low, figure14cseries_high, figure14cseries_mean, figure14cseries_low, figure14dseries_high, figure14dseries_mean, figure14dseries_low, figure14eseries_high, figure14eseries_mean, figure14eseries_low, figure14fseries_high, figure14fseries_mean, figure14fseries_low)
  figure14series <- as.data.frame(figure14series)
  figure14series <- rename(figure14series, quarters = `figure14series.quarter`, figure14ahigh = `figure14aseries_high.V1`, figure14bhigh = `figure14bseries_high.V1`, figure14chigh = `figure14cseries_high.V1`, figure14dhigh = `figure14dseries_high.V1`, figure14ehigh = `figure14eseries_high.V1`, figure14fhigh = `figure14fseries_high.V1`, figure14amean = `figure14aseries_mean.V1`, figure14bmean = `figure14bseries_mean.V1`, figure14cmean = `figure14cseries_mean.V1`, figure14dmean = `figure14dseries_mean.V1`, figure14emean = `figure14eseries_mean.V1`, figure14fmean = `figure14fseries_mean.V1`, figure14alow = `figure14aseries_low.V1`, figure14blow = `figure14bseries_low.V1`, figure14clow = `figure14cseries_low.V1`, figure14dlow = `figure14dseries_low.V1`, figure14elow = `figure14eseries_low.V1`, figure14flow = `figure14fseries_low.V1`)
  
  # Export the series 
  write.csv(figure14series, "figure14series.csv")
}
export_figure14_series()

#---------------------------Figure 13-14 Chart Code-------------------------####
figure13 <- read_csv("[your home directory]/IRF_Figures/OutputData/figure13series.csv") 
figure14 <- read_csv("[your home directory]/IRF_Figures/OutputData/figure14series.csv") 
figure13data <- figure13 %>% filter(quarters<17)
figure14data <- figure14 %>% filter(quarters<17)
theme_if_policystyle <- function() {
  small_margins <- ggplot2::theme(plot.margin = unit(c(2, 2, 2, 3.25), "mm"))
  
  theme_if_policystyle <- ggplot2::theme_classic() +
    theme(
      #text = element_text(family = "Arial"), commented out by m1aem03 due to user reported issues 10/17/2023
      axis.ticks.length = unit(-0.2, "cm"),
      # Remove panel border
      panel.border = ggplot2::element_blank(),
      plot.background = ggplot2::element_rect(fill = "transparent", colour = NA),
      plot.margin = ggplot2::unit(c(0.04, 0.175, 0, 0.175), "in"),
      
      axis.text = ggplot2::element_text(size = 6),
      axis.text.x = ggplot2::element_text(margin = ggplot2::margin(t = 3), color = "black", size = 6), # 10
      axis.text.y = ggplot2::element_text(margin = ggplot2::margin(r = 3), color = "black", size = 6), # 9.5
      axis.text.x.bottom = ggplot2::element_text(size = 6),
      axis.text.y.left = ggplot2::element_text(margin = ggplot2::margin(r = 3), size = 6, hjust = 1), # update 1.3.1 to handle y margin spacing issues
      axis.text.y.right = ggplot2::element_text(margin = ggplot2::margin(l = 3), size = 6, hjust = 1), # 9.5
      axis.title = ggplot2::element_text(size = 6, face = "plain"),
      
      plot.title = ggplot2::element_text(hjust = .5, size = 8, face = "bold"),
      plot.caption = ggplot2::element_text(size = 6, hjust = 0),
      plot.subtitle = ggplot2::element_text(size = 6, hjust = 1, margin = margin(0, 0, 0, 0, "pt")),
      
      axis.ticks = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.line = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.ticks.x = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.line.x = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.line.y = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.ticks.y = ggplot2::element_line(linewidth = 0.3, color = "black"),
      
      legend.text = ggplot2::element_text(size = 8),
      legend.key.size = ggplot2::unit(0.4, 'cm'),
      legend.spacing.x = unit(0.05, 'cm'),
      legend.spacing.y = unit(-0.15, 'cm'),
      legend.key = ggplot2::element_blank(), # removes the legend title
      legend.title = ggplot2::element_blank(),
      legend.background = ggplot2::element_blank()
    )
  
  return(theme_if_policystyle)
}
light_color_set <- c("#afc98d", "#729cc3", "#cec68f", "darkgrey", "orange") 

run_charts <- function(){
#------------------------------Figure13a------------------------------------####
figure13a <- ggplot(data = figure13data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure13amean, ymin = figure13ahigh, ymax = figure13alow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure13amean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(A) Gertler and Karadi (2015)", subtitle = "",) +
  scale_y_continuous(position = "right", breaks = seq(-.02, .05, .01), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.025, .05)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85)) +
  theme(plot.title = ggplot2::element_text(hjust = .5, size = 6, face = "plain"))

figure13a
#------------------------------Figure13b------------------------------------####
figure13b <- ggplot(data = figure13data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure13bmean, ymin = figure13bhigh, ymax = figure13blow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure13bmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(B) Nakamura and Steinsson (2018)", subtitle = "",) +
  scale_y_continuous(position = "right", breaks = seq(-.02, .02, .005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.013, .02)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85)) +
  theme(plot.title = ggplot2::element_text(hjust = .5, size = 6, face = "plain"))

figure13b
#------------------------------Figure14a------------------------------------####
figure14a <- ggplot(data = figure14data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure14amean, ymin = figure14ahigh, ymax = figure14alow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure14amean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(A) Top 0.1%", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.1, .08, .02), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.11, .08)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure14a
#----------------------------Figure14a_ind----------------------------------####
figure14a_ind <- ggplot(data = figure14data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure14amean, ymin = figure14ahigh, ymax = figure14alow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure14amean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(A) Top 0.1%", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.1, .08, .02), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.11, .08)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure14a_ind
#------------------------------Figure14b------------------------------------####
figure14b <- ggplot(data = figure14data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure14bmean, ymin = figure14bhigh, ymax = figure14blow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure14bmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(B) 99-99.9%", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.1, .08, .02), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.11, .08)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure14b
#----------------------------Figure14b_ind----------------------------------####
figure14b_ind <- ggplot(data = figure14data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure14bmean, ymin = figure14bhigh, ymax = figure14blow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure14bmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(B) 99-99.9%", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.1, .08, .02), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.11, .08)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure14b_ind
#------------------------------Figure14c------------------------------------####
figure14c <- ggplot(data = figure14data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure14cmean, ymin = figure14chigh, ymax = figure14clow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure14cmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(C) 90-99%", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.1, .08, .02), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.11, .08)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure14c
#----------------------------Figure14c_ind----------------------------------####
figure14c_ind <- ggplot(data = figure14data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure14cmean, ymin = figure14chigh, ymax = figure14clow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure14cmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(C) 90-99%", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.1, .08, .02), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.11, .08)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure14c_ind
#------------------------------Figure14d------------------------------------####
figure14d <- ggplot(data = figure14data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure14dmean, ymin = figure14dhigh, ymax = figure14dlow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure14dmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(D) 70-90%", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.1, .08, .02), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.11, .08)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure14d
#------------------------------Figure14e------------------------------------####
figure14e <- ggplot(data = figure14data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure14emean, ymin = figure14ehigh, ymax = figure14elow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure14emean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(E) 50-70%", subtitle = "",) +
  scale_y_continuous(position = "right", breaks = seq(-.1, .08, .02), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.11, .08)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure14e
#------------------------------Figure14f------------------------------------####
figure14f <- ggplot(data = figure14data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure14fmean, ymin = figure14fhigh, ymax = figure14flow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure14fmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(F) Bottom 50%", subtitle = "",) +
  scale_y_continuous(position = "right", breaks = seq(-.1, .08, .02), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.11, .08)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.3, 0.85))

figure14f

# Create Figures ####
fig13 <- ggarrange(figure13a, figure13b, ncol = 2, nrow = 1, legend = "bottom", common.legend = TRUE)
fig13a <- ggarrange(figure13a, ncol = 1, nrow = 1, legend = "bottom", common.legend = TRUE)
fig13b <- ggarrange(figure13b, ncol = 1, nrow = 1, legend = "bottom", common.legend = TRUE)
fig14 <- ggarrange(figure14a, figure14b, figure14c, figure14d, figure14e, figure14f, ncol = 3, nrow = 2, legend = "bottom", common.legend = TRUE)
fig14a <- ggarrange(figure14a_ind, ncol = 1, nrow = 1, legend = "bottom", common.legend = TRUE)
fig14b <- ggarrange(figure14b_ind, ncol = 1, nrow = 1, legend = "bottom", common.legend = TRUE)
fig14c <- ggarrange(figure14c_ind, ncol = 1, nrow = 1, legend = "bottom", common.legend = TRUE)
fig14d <- ggarrange(figure14d, ncol = 1, nrow = 1, legend = "bottom", common.legend = TRUE)
fig14e <- ggarrange(figure14e, ncol = 1, nrow = 1, legend = "bottom", common.legend = TRUE)
fig14f <- ggarrange(figure14f, ncol = 1, nrow = 1, legend = "bottom", common.legend = TRUE)

# Save Figures ####
setwd("[your home directory]/IRF_Figures/Figures13_14/")
ggsave(fig13, file = "figure13.eps", height = 3, width = 6, units = "in", device = cairo_ps)
ggsave(fig13a, file = "figure13a.eps", height = 3, width = 3, units = "in", device = cairo_ps)
ggsave(fig13b, file = "figure13b.eps", height = 3, width = 3, units = "in", device = cairo_ps)
ggsave(fig14, file = "figure14.eps", height = 6, width = 9, units = "in", device = cairo_ps)
ggsave(fig14a, file = "figure14a.eps", height = 3, width = 3, units = "in", device = cairo_ps)
ggsave(fig14b, file = "figure14b.eps", height = 3, width = 3, units = "in", device = cairo_ps)
ggsave(fig14c, file = "figure14c.eps", height = 3, width = 3, units = "in", device = cairo_ps)
ggsave(fig14d, file = "figure14d.eps", height = 3, width = 3, units = "in", device = cairo_ps)
ggsave(fig14e, file = "figure14e.eps", height = 3, width = 3, units = "in", device = cairo_ps)
ggsave(fig14f, file = "figure14f.eps", height = 3, width = 3, units = "in", device = cairo_ps)
}
run_charts()